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ABSTRACT 

The three-dimensional, incompressible Navier-Stokes and energy equations 
with the Boussinesq assumption have been directly simulated at a Rayleigh 
number of 3.8 x 10^ and a Prandtl number of 0.76. In the vertical direction, 
wall boundaries were used and in the horizontal, periodic boundary conditions 
were applied. A spectral/finite difference numerical method was used to 
simulate the flow. At these conditions the flow is turbulent, and a 
sufficiently fine mesh was used to capture all relevant flow scales. The 
results of the simulation are compared to experimental data to justify the 
conclusion that the small scale motion was adequately resolved. 
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I. INTRODUCTION 


Direct simulation of turbulent fluid flows is now possible with the large 
vector computers that have become available [1,2], Prediction of low-order 
flow statistics is definitely within current capabilities, and some results 
have already been published which show predictions of small scale turbulent 
features which are consistent with experimental observations [3-5]* The 
current study was undertaken to explore the quality of information that can be 
extracted from a direct flow simulation (DFS) of turbulence on a sufficiently 
fine mesh. 

The turbulent Rayleigh-Benard problem (natural convection) was chosen for 
study since it is a simple turbulent flow for which a good body of 
experimental measurements exists. Moreover, some DFS and large-eddy 
simulations (LES) of this problem have been published albeit on coarser 
meshes. While experimental data do exist, measurements of velocity, where no 
mean flow exists, are difficult. Hence, there is much to be learned about 
turbulent natural convection from an accurate simulation. 

The two requirements for conducting such a study are a high-speed computer 
and an efficient, accurate flow simulation code. The CYBER-205 computer with 
a 16 mega-word memory provides sufficient computation power. This current 
code has been extensively tested, and various versions of it have been used to 
study transition in channel flow [6]. The version used in this study includes 
the addition of the energy equation and a modified vertical momentum equation 
that includes buoyancy consistent with the Boussinesq assumption. 

A simulation of a turbulent flow was then conducted, and these data as 
well as a discussion of the code will be presented in this paper. Though the 
overall goal of this work is an in-depth examination of the quality and type 



2 


of information that can be extracted from such a simulation, the purpose of 
this paper is to document the basic simulation. The simulation results will 
be compared with experimental mean measurements as well as previous DFS 
results. The increased resolution of this work over previous DFS resulted in 
an improvement in the prediction of the Nusselt number; it was sufficiently 
close to experimental results to suggest that in addition to a good prediction 
of the large-scale flow, the small-scale features are accurately represented. 
Grotzbach discusses t\\is connection extensively [7,8]. Comparisons with 
experimental data which are more dependent on the small scale components of 
the flow will also be presented to justify further this conclusion. 


II . RAYLEIGH-BENARD PROBLEM 

The Rayleigh-Benard problem is a simple geometry, laboratory-type problem 
used to study natural convection (Figure 1). Chandrasekhar [9] and Busse [10] 
have described the basic problem and discuss both the stability analysis and 
some experimental results. Krishnamurti [11,12] has summarized much 
experimental data and developed a map showing the qualitative flow at 
different values of Ra and Pr, the principal independent problem parameters 
(defined below). For Pr « 0.76 (air) and Ra = 3.8 x 10 5 the motion is 
turbulent, although it should possibly be qualified as low Reynolds number 
turbulence. Several experimental studies [13-18] and numerical simulations 
[7,19,20] (both LES and DFS) have been completed in the qualitatively similar 
Pr-Ra region. The flow at these values of Pr and Ra consists of a core 
flow (a horizontal layer in the middle 80% of the fluid layer) and a boundary 
region near each plate. The turbulence is statistically homogeneous in the 
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horizontal directions for both layers. In the core the vertical variation of 
most statistical quantities is small. In the boundary layer there is a 
transition from molecular dominated physical processes near the wall to the 
fully turbulent core flow. 

This flow is described by the incompressible Navier-Stokes equations 
modified to include the effect of temperature-induced density variations on 
the buoyancy force (Boussinesq assumption) plus the temperature equation. 
These equations, when non-dimensionalized by a, h, and AT, are 

2 

3u . 3(u,u.) 3P 3 u. 

— - + + p r 1 + Pr Ra TS i3 , (la) 

3t 3Xj 3x^ 3Xj 
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The temperature and pressure in Equations (la,b,c) are the difference between 
the actual temperature, T a , and pressure, P a , and the values due to the static 
temperature gradient only. These are defined as follows: 


T (x,t) = T - x, + T(x,t) 
a — o o — 

= T + T ( x , t ) 
o r — 


and 
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3P a (x» t; ) 




Pr Ra x 3 + 


3P(x, t ) 



The dependent variables in the problem are the velocity components, u^ or 
(u,v,w), the temperature, T, and the fluid pressure, P. The independent 
variables are the spacial coordinates, or _x = (x,y,z), and time, t. The 
indices i = 1 and i = 2 signify the horizontal directions, and i = 3 
denotes the vertical direction* The problem parameters are the thermal 
diffusivity, a; the kinematic viscosity, v; the acceleration of gravity, 
g; the reference fluid density, p q ; the coefficient of thermal expansion, Q; 
a reference temperature (the temperature of the lower plate), T Q ; and the 
temperature difference between the two plates, AT. The dependent and 
independent variables have been non-dimensionalized by a/h (velocity), AT 

n a 

(temperature), p^ a /h (pressure), h (coordinates), and h 2 /ct (time). The 

3 

Rayleigh number, Ra = — , an( j Prandlt number, Pr = v/a, are the 
principal non-dimensional problem parameters. 


III. NUMERICAL METHOD 

The Fourier finite difference algorithm developed by Moin and Kim [21] 
for their large-eddy simulations of turbulent channel flow has been applied, 
with modifications, to the present direct simulation of turbulent Rayleigh— 
Benard flow. This Rayleigh-Benard Fourier-finite difference Method (RBFFDM) 
is an unsplit method on a grid staggered (for the pressure variable and the 
continuity equation) in the vertical (z) direction. Fourier collocation is 
used for the spacial discretization in the x and y directions whereas in 
the vertical direction second-order finite-differences are employed on the 
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non-uniform grid, 


z = (l - cos(kTr/N ))/ 2, 


k = 0,1, ••^N . 
’ * z 


The time discretization is Crank-Nicolson for the viscous and conductive 
terras and backward Euler for the pressure gradient term. The advection and 
buoyancy terms are handled by a third-order Adams-Bashf orth method. 

The implicit part of the algorithm requires, for each pair of horizontal 
Fourier wavenumbers, the solution of 2 real, block tridiagonal systems 
(involving the velocities and pressure) and, independently, 2 real, scalar 
tridiagonal systems (for the temperature). The block-tridiagonal equations 
were scaled as described by Zang & Hussaini [6] for their Fourier-Chebyshev 
version of the corresponding channel flow algorithm. Pivoting has proven to 
be unnecessary for this system. The block-tridiagonal solution algorithm 
takes advantage of the many zero elements which occur. Vectorization of this 
phase of the algorithm is achieved by solving for many pairs of Fourier 
wavenumbers at the same time. 

Equations (la,b,c) were solved on the region, 0<x<A,0<y<A, 

x y 

0 < z < 1. Under the present scaling of the vertical direction, the lengths 
A x and A y correspond to the aspect ratios of the two horizontal directions 
to the vertical one. The boundary conditions at the lower and upper walls, 
z - 0 and z = 1, are the conventional no-slip and no temperature jump 
conditions. In the horizontal directions, periodic boundary conditions are 
assumed. These aspect ratios must be large for reasonable correspondence with 
experiments. The aspect ratios are related to the resolution by A^ = Ax 
and Ay = Ay. The computer memory limitation places an upper bound on 
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N = N x Ny N z ; therefore, to have a large A x ^ along with sufficiently small 
Ax i? some compromise is necessary. The values chosen, A x = 4 and Ay = 2, 
allow the available computer memory to be used for better small scale 
resolution. This will be further discussed in Section IV. 


The RBFFDM code has been implemented on a CDC Cyber 205 with 2 pipes and 
16 million 64-bit words of main memory. For each grid point 13 variables were 
stored. Additional storage equivalent to 7 variables per grid point was used 
to facilitate vectorization. A total of 11 million words was used for the 


simulation on the 128 x 64 x 64 grid. Vector lengths for the explicit portion 
of the algorithm were between 4 and 40 thousand. Typically, one-fourth of the 
implicit equations were solved together. The vector lengths here were roughly 
1,000. The linked triad feature was heavily exploited. A single time-step 
required 6.8 seconds of CPU time and no 1/0 time since the job was run 
entirely within the central memory. The sustained speed of the calculation 
was 100 MFLOPS. 


The start-up phase of the calculation took 2100 time-steps and the data 
collection an additional 5600 steps. A total of 12 hours of CPU time was 


required for the data collection. This includes the time for some preliminary 
diagnostics. The Courant number, defined as the maximum over the grid of the 


quantity 


M + M + M 

Ax Ay Az 


At , 


ranged between 0.19 and 0.26 and averaged 0.23 for the data collection phase. 
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IV. COMPARISON WITH PREVIOUS WORK 

The turbulent flow of the Rayleigh-Benard problem is assumed to be 
homogeneous in both horizontal directions as well as statistically steady in 
time (after a start-up period). Experimental data are usually presented using 
some combination of a long-time average as well as a spacial average in one or 
both horizontal directions. The simulation results presented as horizontal 
averages, < >, have been averaged over both the x and y directions (except 
for the 1-D x-spectra which were averaged in y only). In addition, they 
have been time averaged over a time period equal to 10/W C , where 
W c = (Nu Pr Ra)*/-^. (See Deardorff [22] for a discussion of W c , the scaling 
velocity for the large eddies.) This period, which should consist of several 
large eddy turn-over times, was found adequate by Eidson [19]. Volume 
averages were time averaged as well. All the simulation results presented 
below are both horizontally and temporally averaged unless otherwise 
specified. 

Prediction of the average vertical heat flux, the Nusselt number in non- 
dimensional form, is an important result of any natural convection study. 
Previous simulations have predicted values of the Nusselt number, Nu, which 
are slightly higher than those measured experimentally. Grotzbach [8] has 
discussed extensively this discrepancy and has shown that inadequate 
resolution is partially responsible. In Table 1, the results from both 
simulations and experiments are shown. The prediction of the present study 
lies at the the upper range of experimental measurements and below the Nusselt 
numbers predicted by previous, coarser-grid simulations. In the present 
simulation, Nu was calculated at each z level. The average value is 
reported in Table I. The variation with z was small (approximately ± 0.1) 
except very near the lower wall where Nu increased to 7.0. 
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The aspect ratio of the horizontal to vertical boundary lengths also is 
known to affect Nu [8,15]. The values of A x = 4 and Ay = 2 are smaller 
than the values of 4 to 7 suggested by experimentalists as the minimum for 
removing significant side boundary effects. Although the aspect ratio effect 
is not negligible, especially on Nusselt number [8], a large aspect ratio was 
foregone in favor of better small scale resolution in view of the goal of this 
study: to resolve eddies down to nearly the dissipation scales. 

The relative temperature, T r , vertical profile is compared to data of 

8T 

studies by Chu and Goldstein [16] in Figure 2. A line with slope, — — — = Nu, 
is drawn in this figure. From this one can estimate the conductive layer 
thickness, 6^. This will be used later in a more extensive examination of 
the temperature data near the wall. 

In Figures 3 and 4, the vertical dependence of the velocity and 
temperature RMS values are compared to the simulation of Grotzbach [7] and the 
measurements of Deardorff and Willis [13] (slightly larger Ra). In Table 1 a 
comparison for the centerline values is presented for a wide range of studies. 
The present improved resolution reduces slightly the predicted levels of the 
RMS temperature over previous simulations. However, no systematic change in 
the RMS velocity levels with the improved spacial resolution was observed. 
Both trends are consistent with Grotzbach's results. Considering the 
variation In experimental values and the uncertainty in aspect ratio effects, 
the results appear quite satisfactory. 

The wT correlation coefficient, C(wT), is constant for 0.2 < z < 0.8, 
giving a value of 0.71 (Figure 5). Both Grotzbach [7] and Eidson [20] 
previously obtained a value of approximately 0.67. Deardorff and Willis 
measured approximately 0.60 for a slightly higher Ra, but C(wT) should 
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decrease with increasing Ra. Near the wall all four studies differ. 
Deardorff and Willis measured a significant drop as the wall was approached, 
but warned that since the numerator and denominator of C(wT) become small, 
their results are uncertain. The increase in C(wT) near the wall found in 
the current study was not observed in either previous simulation, but neither 
of these had sufficient resolution in this region. Note that only the data of 
GrotzbaclTs case 9 with 16 vertical grid points was available for 
comparison. Also, horizontal averages of the uv, uw and vw correlations 
were calculated. They were all approximately zero as would be expected for 
turbulence homogeneous in the horizontal directions. 

In Figure 6, the several terms in the horizontally averaged kinetic 
energy equation are plotted versus z and compared with experimental data at 
a slightly larger value of Ra. These terms are 


3<E> 

IT" 


= Pr Ra <wT > - 
a 

Production 


3 

"ST 


<w(E + P )> 
a 

Diffusion 


3u 3u. 

" Pr 

J J 


Pr 


3 2 <E> 

3 ? 


Dissipation Molecular Transfer 


where 

E = 7 u i u i * 

Since the production term in the core is equal to (Pr Ra Nu), the 
experimental data of Deardorff & Willis [13] at Ra = 6.3 x IQ"* and Nu = 5.7 


would be expected to be about 40% higher than the current simulation. For 



10 


clarity, only a few key values from Deardorff and Willis are shown to 
demonstrate the general agreement between experiment and simulation. The four 
terms are all essentially constant and approximately 40% below the 
experimental data for the core region. In the boundary region the variation 
with z is similar for both experiment and simulation. The diffusion term 
(for the experiment and simulation) peaks at a z value of 0.035, and the 
molecular diffusion term of the simulation peaks at a slightly higher z. The 
experimental point plotted at a low z value is at these peaks for these two 
terms. Actually, one would expect the peaks for the lower Ra to be at a 
slightly larger z since the boundary region thins with increasing Ra. From 
a closer inspection of the data, the molecular dissipation also can be seen to 
have a small positive value between z = 0.1 and z = 0.2 as was found by 
Deardorff and Willis. The production and dissipation terms have the same 
general shape in both studies. The volume averages of the production and 
dissipation terms were 1.641 x 10 6 and 1.637 x 10 6 , respectively. Also, the 
volume average of the molecular transfer term was 9.0 x 10^. 

One advantage of the simulation technique is that quantities which are 
experimentally difficult to measure can be easily calculated. Figures 7 to 9 
give three examples. In Figure 7 the absolute value of each of the three 

components of vorticity is shown. As expected, the x and y components are 

large near the wall due to the creation there of a boundary layer by the large 

eddies. Near the center the flow is more isotropic. The ratio of the volume 
average of the horizontal to vertical vorticity was 5.0 and 4.3 for the x 

and y directions, respectively. This is higher by a factor of 2 than in the 
previous DFS by Eidson [19]. Figure 8 shows that the only significant 
velocity skewness is for the w component near the wall. A negative value 
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for the w component is reasonable since the fluid particles with negative 
velocity (near the lower wall) come from the core where there is more velocity 
variation due to the turbulent cascading process. The particles with a 
positive velocity originate near the wall where the motion is damped and they 
have a more uniform value. Recall that the horizontal average of the w 
velocity equals zero. The flatness profile of the velocities is shown in 
Figure 9. For the horizontal components, these approach 3 in the core, a 

value which is similar to that in channel flow turbulence away from the 
wall [21]. The larger flatness factor near the wall suggests that the flow is 
more intermittent in this region. 

The wT cross-spectra for several z values are plotted in Figures 10 
to 13. These are 1-D spectra in the x-direction. They have been normalized 
so that they integrate to <wT>. The wavenumber, k, is defined such that 2ttkx 
is the argument of the periodic expansion functions. Notice that the spectra 
decrease with k except near the wall where they are roughly constant with 
k. This trend was found also in the u, w and T spectra but not in v. 
Near the center the most energy would be expected in the largest available 
scales for the small values of A in this study. Near the wall, smaller 
horizontal structures would be expected, and thus significant energy at a 
wavenumber greater than the minimum is a reasonable result. Since this is low 
Ra turbulence, an equilibrium region would not be expected, and indeed none 
was found here. The drop-off with wavenumber is rapid due to the importance 
of viscous effects at all the scales. 

To evaluate the simulation further, a more sensitive measure of the 
small scale features of the flow field, especially in the wall region, was 
sought. Carroll [14] has made direct instantaneous measurements of the 
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temperature and the vertical temperature gradient and calculated the mean, 
standard deviation, and skewness. Carroll's data are compared with the 
results of this simulation in Figures 14 to 20. The data of Carroll in these 
figures are a mean curve (drawn by Carroll) through the data points. The bars 
show the range of data scatter. The mean of the relative temperature 

(simulation data shown in Figure 2) and its vertical derivative (Figure 18) 
show good agreement with Carroll's data. The higher order statistics give a 
variation with z which is similar to the experimental results although the 
magnitudes differ in some cases. In Table 1, a significant variation in 

experimental measurements of the RMS of T r> R(T r ) , between different 

experimental studies is apparent. Carroll used a stationary probe, and this 
may have reduced the fluctuations coming from the direction of the probe and 
lowered the statistical levels. Also Carroll measured the same RMS 

temperature level for Nu ranging from 5.5 to 14.0. Both experiments [13,15] 
and simulations [20] have shown that R(T r ) increases with decreasing Nu 

(or Ra). Since the magnitude of the data is in question in some cases, the 
conclusions drawn from the simulation/experiment comparison are based mainly 
on changes of slope with vertical distance from the wall. These variations 
correspond to differences in the dominant physics in the different layers, and 

the simulation was able to predict these layers in agreement with the 

experimental results. 

A four-layer (or three-layer with the middle divided into 2 sub-layers) 
model is consistent with both the experimental and simulation* results . These 


are: 
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(i) a conduction dominated layer 


0 < Nu*z < Sj* 
0.1 < 5* < 0.3 


(ii) an inner boundary layer 


(iii) an outer boundary layer 


(iv) the core 


6 *< Nu*z < 6 2 * 

<$ 2 * - 0.5 

6 2 * < Nu*z < 6 3 * 

* 

6 3 - 1.5 

6 3 * < Nu*z < Nu/2 . 


These layers are all well resolved by the simulation as they contain 8, 4, 9, 

£ 

and 11 grid points respectively (assuming 6^ = 0.2). 

A 

Determination of 6^ is uncertain. As was mentioned previously (Figure 

2), divergence of the temperature profile from a straight line with a 

normalized slope equal to Nu (the temperature gradient at the wall) gives 
* * 

6, =6 = 0.26. Another measure is that the wT correlation rises from 1 % 

l c 

to 10% of its core value near Nu*z » 0.1. 

The inner boundary layer is characterized by the transition from 

molecular to turbulent dominated processes. Both experiment and simulation 
3T r 

show that has been reduced by 50% between Nu*z of 0.4 and 0.5 (Figure 

18). In the inner boundary layer, R(T r ) reaches a maximum (Figure 14), the 
skewness factor for T r , S(T r ), changes sign (Figure 16), and the flatness 
factor for T r , F(T r ), reaches a minimum (Figure 17). The skewness, S(T r ), 
continues to increase reaching a positive maximum in the outer boundary layer. 
The above trends, shown in both studies, suggest significant changes in the 
flow character in these layers. 
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In the outer layer, the flow transitions to a core region where the 

large eddies (the most efficient heat transfer agents) carry the heated fluid 

3T 

from the bottom to the top. At Nu*z w 1.5, the temperature gradient, — — — t 
has been reduced to 1% of its maximum value (Figure 18). The skewness, S(T r ), 
decreases in the outer layer and approaches zero near the center of the core 
(Figure 16). Also note, in the outer layer R(T r ) varies as z~^^ compared 
with z +1 in the inner layer, a trend which is found in both studies (Figure 
15). The simulation predicts a transition to a constant level for R(T r ) in 
the center in agreement with Deardorff and Willis [13] (Figure 4). However, 
Carroll's data suggest that R(T r ) varies as z“*/3 i n core as we ll as 

the outer layer. 

A possible discrepancy between the simulation results and Carroll's data 
3T 

is seen in the R( ^ z ) plot (Figure 19). At values of Nu*z between 0.4 
and 1.0 the simulation shows a hump. However, Carroll's choice of an "average 
curve" has significant scatter only on top of the curve. Re-drawing the data 
fit through this scatter would result in a similar hump in the outer boundary 
layer region. Another example, where the simulation results suggest that the 
"averaged curve" through sparse, scattered data is not correct, is Carroll's 
choice for the temperature skewness profile for Nu*z <0.1 (Figure 16). He 
suggests that the curve should continue decreasing to S(T r ) = -1.5 and then 
rise to zero at Nu*z *• 0.01. However, only a few data points which gave 
significant scatter were measured in this region. From an examination of the 
data in Carroll's paper, it is not inconsistent to draw the experimental data 
fit to turn between Nu*z = 0.1 and 0.2 and be parallel to the simulation 
results somewhere between the dashed lines shown in Figure 16. One 
discrepancy that cannot be explained by data scatter is the comparison between 
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studies of sf-^ — 1 away from the wall. The simulation predicts that sf-* — ] 
returns to zero in the outer layer and remains zero in the core (Figure 20). 
However, the experiments show a return to negative skewness in the core. 

The final comparison is for the temperature gradient versus Nu*z near 
the wall (Figure 18). Various power laws have been hypothesized which predict 
a linear slope for the data plotted in log-log form. Carroll [14], Businger 
[23], and Monin and Yaglom [24] discuss the various theories. Carroll's data 
show that only for the range 0.4 < Nu*z <1.0 does a simple power law seem 
reasonable and even then the data have a slow change in slope in this 
region. The current simulation fits Carroll's data nicely and thus also 
predicts a slope of -2 in this range of z. The value of -2 is predicted by 
the theory of Malkus [25]. 

From these comparisons, it can be seen that the simulation data 
represent a turbulent realization which has a good agreement with 
experimental, averaged measurements. The agreement, at least in prediction of 
sub-layers using the variance and skewness data, shows that the small scale 
turbulent features are accurately represented near the wall. This is 
important because a complete understanding of the Rayleigh-Benard convection 
requires studying the thin boundary layer which exists as a result of the 
large scale motion in the core. From isothermal, flat plate boundary layer 
studies, it is known that significant small scale events occur in this type of 
region which are important to the global flow. 

Some of the features that we have observed in this calculation were also 
found in a reference pointed out to us by a referee [26]. However, we have 
not yet performed a detailed comparison with this data. 
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V. CONCLUDING REMARKS 

From the comparison with the gross properties of the experimental data, 
it can be concluded that the completed computer simulation is a reasonably 
accurate realization of a turbulent flow even down to the dissipation scales. 
Moreover, the detailed database has enabled us to identify a feature in the 
skewness of the vertical temperature derivative which has been overlooked in 
the experiment due to data scatter. 

Obviously, a far more detailed analysis of the simulation results is in 
order. This simulation is of sufficient accuracy to study the small scale 
structure of a natural convection flow in a similar manner to the work 
previously done for channel flows [4,21]. Moreover, it can be used to help 
interpret experimental measurements which are typically limited to a few 
spacial points and fewer than all three velocity components. Finally, we feel 
this database can be used to evaluate turbulent models and other theoretical 
ideas. 
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Table 1. Horizontal Averages at the Vertical Center 


Direct Simulation 

Grid 

Ra 

Nu 

U RMS 

W RMS 

T 

r RMS 

Current Study 

128x64x64 

3.8xl0 5 

6.6 

52 

116 

0.08 

Grotzbach* 


3.8xl0 5 





(case 7) 

16x16x16 


7.8 

— 

109 

— 

(case 9) 

32x16x16 


7.4 

50 

107 

0.08 

(case 14) 

64x32x32 


6.9 

— 

109 

0.08 

Eidson 

64x64x16 

3.8xl0 5 

8.1 

70 

112 

0.11 


Experiment 






Deardorff & Willis 

6.8xl0 5 

5.8 

65 

120 

0.08 

Carroll 

3.8x10^ 

5.9* 

— 

— 

0.05* 

Fitz jarrold 

3.8xl0 5 

6.2* 

80* 

120* 

0.05* 

Chu & Goldstein 

3.8xl0 5 

6.5* 

— 

— 

— 

Goldstein & Chu 

3.8xl0 5 

5.4* 

— 

— 

— 


Calculated from curve fit of data over a range of Ra 
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Figure 3. Comparison with experiments and previous 

simulations of the RMS of the velocity components 
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Figure 5. wT correlation versus height. 
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Figure 6. Terms in the kinetic energy equation versus height. 




Figure 7. RMS of the three vorticity components versus height. 
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Figure 8. Skewness factor of the three velocity components versus height 
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Figure 9. Flatness factor of the three velocity components versus height 
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Figure 16. Comparison with experiment of the skewness of the 
relative temperature versus normalized height. 
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Figure 18. Comparison with experiment of the normalized relative 
temperature derivative versus normal ized height. The 
slope shown is from the theory of Malkus. 
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Figure 20. Comparison with experiment of the skewness factor 

for the normalized relative temperature derivative. 
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